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Abstract 

After sketching the basic principles of renewal theory we discuss the 
classical Poisson process and offer two other processes, namely the renewal 
process of Mittag-Leffler type and the renewal process of Wright type, so 
named by us because special functions of Mittag-Leffler and of Wright type 
appear in the definition of the relevant waiting times. We compare these 
three processes with each other, furthermore consider corresponding renewal 
processes with reward and numerically their long-time behaviour. 
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1. Introduction 

It is well known that the Poisson process (with and without reward) plays a 
fundamental role in renewal theory. In this paper, by means of functions of 
Mittag-Leffler and Wright type we provide a generalization of and a variant 



to this classical process and construct interesting subordinated stochastic 
processes of fractional diffusion. 

The plan of the paper is as follows. In Section 2 we recall the basic 
renewal theory including its fundamental concepts like waiting time between 
events, the survival probability, the renewal function. 

If the waiting time is exponentially distributed we have the classical 
Poisson process, which is Markovian: this is the topic of Section 3. 
However, other waiting time distributions are also relevant in applications, 
in particular such ones with a fat tail caused by a power law decay of their 
density. In this context we analyze, respectively in Sections 4 and 5, two 
non-Markovian renewal processes with waiting time distributions described 
by functions of Mittag-Lefner and Wright type, that exhibit a similar power 
law decay. They both depend on a parameter (3 6 (0, 1) related to the 
common exponent in the power law. In the limit (3=1 the first becomes 
the Poisson process whereas the second goes over into the deterministic 
process producing its events at equidistant instants of time. 

In Section 6, after sketching the basic differences between the renewal 
processes of Mittag-Lefner and Wright type, we compare numerically their 
survival functions and their probability densities in the special case (3 = 1/2 
with respect to the corresponding functions of the classical Poisson process. 

In Section 7 we discuss the general renewal process with reward 
(also called the compound renewal process), interpret it physically as a 
continuous time random walk (CTRW), and review how the resulting sojourn 
probability density function, evolving in time according to an integral 
equation, can analytically be represented as an infinite series. We then 
consider the compound renewal processes of Mittag-Lefner and Wright type 
including the limiting case (3 = 1. In Section 8, we recall the conditions 
of the well-scaled transition to the diffusion limit according to which the 
CTRW integral equation of our compound processes reduces to the time 
fractional diffusion equation (TFDE). For all three cases (Poisson, Mittag- 
Lefner, Wright), by numerically summing the corresponding series for the 
cumulative probability function, we produce plots of its behaviour in time. 
These plots nicely illustrate that with increasing time the behaviour of time- 
fractional diffusion is approached. Finally, concluding remarks are given in 
Section 9. 
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2. Essentials of renewal theory 



For the sake of reader's convenience, we present a brief introduction to the 
renewal theory by using our notation. For more details see e.g. the classical 
treatises by Cox [8], Feller [llj . and the more recent book by Ross [41] . 
We begin to recall that a stochastic process {N(t), t > 0} is said to be a 
counting process if N(t) represents the total number of "events" that have 
occurred up to time t. The concept of renewal process has been developed 
as a stochastic model for describing the class of counting processes for which 
the times between successive events, Ti,T 2 , . . ., are independent identically 
distributed (iid) non-negative random variables, obeying a given probability 
law. We call these times waiting times (or inter-arrival times) and the times 

k 

to = 0, t k = J2 T i, k > 1 - (2- 1 ) 

3=1 

renewal times. That is, to = is the starting time, t\ = T\ is the time of 
the first renewal, t% = T\ + T 2 is the time of the second renewal, and so on; 
in general tk denotes the time of the feth renewal. 
Let the waiting times be distributed like T and let 

*(t) := P(T < t) (2.2) 

be the common probability distribution function. Here P stands for 
probability. We assume <3?(i) to be absolutely continuous, so that we can 
define its probability density function (p(t) as 

0(t) = ^*(t), *(t) = J\{t')dt' . (2.3) 

We recall that <j){t) > with / °° <j>(t) dt = 1 and is a non-decreasing 
function in R + with $(0) = 0, $(+00) = 1. 

Let us remark that, as it is popular in Physics, we use the word 
density also for generalized functions that can be interpreted as probability 
measures. In these cases the function <$(i) may lose its absolute continuity. 
Often, especially in Physics, the probability density function is abbreviated 
by pdf, so that, in order to avoid confusion, the probability distribution 
function is called the probability cumulative function (being the integral of 
the density) and abbreviated by pcf. When the nonnegative random variable 
represents the lifetime of technical systems, it is common to call <£(i) the 
failure probability and 

/oo 
<j)(t') dt' = 1- $(i) , (2.4) 
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the survival probability, because and are the respective probabil- 
ities that the system does or does not fail in (0,t\. These terms, however, 
are commonly adopted for any renewal process. 

As a matter of fact the renewal process is defined by the counting process 



N(t) is thus the random number of renewals occurring in (0, t]. We easily 
recognize that = P(N(t) = 0) . For an example of a renewal process, 
suppose that we have an infinite supply of light-bulbs whose lifetimes are 
iid random variables. Suppose also that we use a single light-bulb at a time, 
and when it fails we immediately replace it with a new one. Under these 
conditions, {N(t), t > 0} is a renewal process when N(t) represents the 
number of light-bulbs that have failed by time t. 

Continuing in the general theory, we set i*i(t) = 3>(i), fi(t) = 4>(t), and 
in general 



F k (t) := P (t fc = Ti + . . . + T k < t) , f k (t) = -F k (t) , k > 1 . (2.6) 



F k (t) is the probability that the sum of the first k waiting times does not 
exceed t, and f k (t) is the corresponding density. F k (t) is normalized because 



k random variables each of which is finite with probability 1 is finite with 
probability 1 itself. 

We set for consistency Fo(t) = &(t), the Heaviside unit step function 
(with 6(0) := 6(0+)) so that F (t) = 1 for t > 0, and /„(t) = S(t), the 
Dirac delta generalized function. 

A relevant quantity related to the counting process N(t) is the function 
v k (t) that represents the probability that k events occur in the closed interval 
[0,t]. Using the basic assumption that the waiting times are i.i.d. random 
variables, we get, for any k > 0, 



We note that for k = we recover vo(t) = ^(t). 

Another relevant quantity is the renewal function m(t) defined as the 
expected value of the process N(i), that is 




for < t < h 
for t > h . 



(2.5) 




v k (t) := P (N(t) = k) = P(t k <t, t k+1 >t)= ( f k (t') V(t - 1') dt' . (2.7) 



Jo 



rn 



(t) :=E(N(t)) = (N(t)) = Y,P(tk<t). 



(2.8) 



fc=i 
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Thus this function represents the average number of events in the interval 
(0,t] and can be shown to uniquely determine the renewal process [41] . It is 
related to the waiting time distribution by the so-called Renewal Equation, 



m(t) = $(£)+/ m(t-t')(j){t')dt' = [1 + m{t - t')] cj)(t') alt' . (2.9) 
Jo Jo 

If the mean waiting time (the first moment) is finite, namely 

POO 

p:=(T)= t <t>(t) alt < oo , (2.10) 



it is known that, with probability 1, t^/k — > p as k — > oo , and N(t)/t — > 1/ p 
as t —* oo . These facts imply the elementary renewal theorem, 

m(t) 1 

— ^ -> - as i^oo. 2.11) 

We shall also consider renewal processes in which the mean waiting time 
is infinite because the waiting time laws have fat tails with power law 
asymptotics: 

Ht) ~ , *(*) ~ jjp, fort -do, 0</5<l, Aoo >0. (2.12) 

It is convenient to use the common * notation for the Laplace 
convolution of two causal well-behaved (generalized) functions f(t) and g(t), 

f{t') g(t - t') alt' = (f*g) (t) = (g * /) (t) = f f(t - t') g(t') alt' . 



Being /& (i) the palf of the sum of the k Hal random variables T\ , . . . , T& with 
pcif (/>(£) , we recognize that fk(t) is the A;- fold convolution of <p{t) with itself, 

= (/*)(*), (2-13) 

so that Eq. (2.7) simply reads: 

w fc (t) := P (JV(t) = fc) = (* * <^ fc ) (t) . (2.14) 

We note that in the convolution notation the renewal equation (2.9) reads 

m(t) = *(*) + (m * (f)){t). (2.15) 
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The presence of Laplace convolutions allows us to treat a renewal process 
by the Laplace transform. 

Throughout we will denote by f(s) the Laplace transform of a sufficiently 
well-behaved (generalized) function f{t) according to 

„ /-+oo 

£{f(t);s} = f(s)= / e- st f(t)dt, s>s , 
Jo 

and for 5(t) consistently we will have S(s) = 1 . Note that for our purposes 
we agree to take s real. 

We recognize that (2.13)-(2.14) reads in the Laplace domain 



fk(s) = \<f>(8)] k , V k (s) = \<f>(s)] k *(S) , 



where, using (2.4), 



1 - <f>(s) 



Then, in the Laplace domain, the renewal equation reads 



m(s) = <&(s) + m(s) (j)(s) , with $(s) 



from which 



m(s) 



1 - cf>(s) 



sm(s) 



1 + s m(s) 



(2.16) 
(2.17) 

(2.18) 
(2.19) 



3. The Poisson process as a renewal process 

The most celebrated renewal process is the Poisson process (with parameter 
A > 0). It is characterized by a survival function of exponential type: 



(t) = e 



-At 



t > 0. 



(3.1) 



As a consequence, the corresponding density for the waiting times is 
exponential as well: 

0(t) = \e~ xt , t>0, (3.2) 



with moments 



(3.3) 
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and the renewal function is linear: 

m(t) = Xt, t>0. (3.4) 

The Poisson process is Markovian because the exponential distribution is 
characteristic for a process without memory. We know that the probability 
that k events occur in an interval of length t is the celebrated Poisson 
distribution with parameter At, 

v k (t):=P(N(t) = k) = ^e~ xt , t>0, k = 0,1,2, .. . . (3.5) 

The probability distribution related to the sum of k iid exponential random 
variables is known to be the so-called Erlang distribution (of order k > 1). 
The corresponding density (the Erlang pdf) is thus 

^) = A {/r^ e ~ A< > ^ = 1 ,2,..., (3.6) 

and the corresponding Erlang pcf is 

F k (t) = f f k {t>) dt' = l- J2 ^ e- A * = £ Wf- e~* , t > . (3.7) 

In the limiting case k = we recover fo(t) = 5(t) and Fo(i) = @(t). 

The results (3.4)-(3.7) can be easily obtained by the Laplace transform 
technique. In fact, starting from the Laplace transforms of the probability 
laws (3.1)-(3-2) 

*W = r J r , ^) = i^-. ( 3 - 8 ) 

A + s A + s 

and, using (2.16)-(2.19), we have, 

Ms) = r n = 4 . Ms) = ${')]" = A l +1 , (3-9) 

s 1 — <p(s) 
hence (3.4)-(3.5), and 



s 2 ' KV 7 Lrv 7J v ; (A + S ) fe +! 



hence (3.6)-(3.7). 
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4. The renewal process of Mittag-Leffler type 

A "fractional" generalization of the renewal Poisson process has been 
recently proposed by Mainardi, Gorenflo and Scalas 29J. Noting that the 
survival probability for the Poisson renewal process (with parameter A > 0) 
obeys the ordinary differential equation (of relaxation type) 

^¥(t) = -A* (t) , t>0; *(0+) = l. (4.1) 

the required generalization is obtained by replacing in (4.1) the first 
derivative by the fractional derivative (in Caputo's sense] of order (5 E (0, 1]. 
We thus write, taking for simplicity A = 1, 

t Z>£tt(i) = -¥(t), t>0, < /3 < 1 ; *(0+) = l. (4.2) 

The solution of Eq. (4.2) can be obtained by using the technique of the 
Laplace transforms. We have for t > : 

*(t) = Ep{r-lP) , from #( S ) = ^-A, 0</3<l, (4.3) 

hence : 

<t>{t) = ~^*(t) = -^jE^-tP), corresponding to 0(a) = j^p- (4.4) 

Here Ep denotes the Mittag-Leffler functional of order (3. We call this process 
the renewal process of Mittag-Leffler type (of order ($) . 

The Caputo derivative of order j3 £ (0, 1] of a well-behaved function f(t) in R + is 

1 f W (r) 



tD* f(t) := 



V{l-0)J o (t-rY 



dr, O<0<1, 



f t f(t), P = l. 



Its Laplace transform turns out as 

C{ t D^f(t)-s} = s"f(s)-s''~ 1 f(0 + ). 

For more information on the theory and the applications of the Caputo fractional derivative 
(of any order (3 > 0), see e.g. ® El EH SS3] ■ 

2 The Mittag-Leffler function with parameter (5 is defined as 

°° 

It is an entire function of order j3 and reduces for (3 = 1 to exp(z) . For detailed information 
on the functions of Mittag-Leffler type the reader may consult e.g. [TOl [TBI l24l [28l l40l [43] 
and references therein. 
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Hereafter, we find it convenient to summarize the most relevant features 
of the functions ^f(t) and 4>{t) when < (3 < 1 . We begin to quote their 
expansions in power series of t& (convergent for t > 0) and their asymptotic 
representations for t — ► oo, 

ifro T(/?n + l) it tP ' K ' 

MA 1 IP ^ S m(/?vr) r(/3 + l) 

We recognize that for < (3 < 1 both functions VP(£), keep the complete 
monotonicitjjf], characteristic for the Poissonian case f3 = 1. 

In contrast to the Poissonian case, in the case < (3 < 1 the mean waiting 
time is infinite because the waiting time laws no longer decay exponentially 
but exhibit power law asymptotics according to (2.12) where the constant 
Aoo derived from (4.5)-(4.6) is 

A OQ = T(P+l) sin(/^)/7r = 13 T{(3) sin(/?vr)/7r . (4.7) 

As a consequence the process turns out to be no longer Markovian but of 
long-memory type. 

The renewal function of this process can be deduced from the Laplace 
transforms in (2.19) and (4.4); we find 

1 iP 
m(s) = ^, hence m{t) = — — , t>0, < < 1 . (4.8) 

Thus, for (3 < 1 the renewal function turns out super-linear for small t and 
sub-linear for large t. 



Complete monotonicity of a function f(t) means, for n = 0,1,2,..., and t > 0, 

d n 

(—l) n -—f(t) > 0, or equivalently, its representability as (real) Laplace transform of a 
non-negative function or measure. Recalling the theory of the Mittag-Lemer functions of 
order less than 1, we obtain for < f3 < 1 the following representations, see e.g. |16| . 

(t) ~ TV J r 2 ' 3 + 2r' 3 cos(/37r) + l ' *~ ' 
„, Bin (/for) f°° r e~ rt 

4>{t) = ii— i / —. —, —— dr, t > . 

^ w n J r 2 ' 3 + 2r' 3 cos(/3tt) + 1 
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10 2 10 1 10° io' io 2 "io 2 10~' 10° io' io 2 

t t 

Figure 1: The functions (left) and <j>{t) (right) versus t (1CT 2 < t < 10 2 ) 
for the renewal processes of Mittag-Lemer type with (3 = 0.25,0.50,0.75, 1. 

For the generalization of Eqs (3.5)-(3.7), concerning the Poisson and the 
Erlang distributions, we give the Laplace transform formula 

£{tP k Ej*\-tP); 8 }= (1 ^fc + i > ^ >0 ' ^ = 0,1,2,.... (4.9) 

with E£ >(z) := —^E p (z) , that can be deduced from the book by Podlubny, 
see (1.80) in iOl. Then we get, with < < 1 and k = 0, 1, 2, . . ., 

v k (t) = ^Ef\- t P), from v k (s) = 9(s) [^)] fc = ^ , (4.10) 

as generalization of the Poisson distribution (with parameter t), what we 
call the (5-fractional Poisson distribution. Similarly, with < (3 < 1 and 
k = 1, 2, . . ., we obtain the (3-fractional Erlang pdf (of order k > 1): 

/fe(*)=^77^w4 fromfe( S ) = [^)f= (4.11) 

and the corresponding (3-fractional Erlang pcf: 

F k (t) = / /*(*') dt' = 1 - J2 -r Ef\-t?) = J2— Ef\-t?) . (4.12) 
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5. The renewal process of Wright type 



A possible choice for obtaining an analytically treatable variant to the 
Poisson process has been suggested by Mainardi et al. [33]. It is based on 
the assumption that the waiting-time pdf <p(t) is the density of an extremal, 
unilateral, Levy stable distribution with index (3 S (0, 1), which exhibits, as 
we shall show, the same power law asymptotics as the corresponding pdf of 
the previous renewal process. In this case, however, the transition to the 
limit /3 = 1 is singular, and the Poisson process is no longer obtained. Now 
we have for t > 0, 

i£(t) = { 1 " ""ft 1 ' </ 3 < 1 ' from V(s) = 1 ~ 6 S , (5.1) 

1 e(t-i), (3 = i, s 




~^)' °<P< 1 > fromS(s) = e~ fi/3 , (5.2) 
P = l, 

where \ and &-/3,o denote Wright functional). 

In view of the presence of the Wright function in (5-l)-(5.2), we call this 
process the renewal process of Wright type. 

Hereafter, taking < (3 < 1, we quote for ^(t) and 4>{t) their expansions 
in powers series of t~^ {convergent for t > 0) 

•O^EH)-'^^, </»<!, (5.3) 

n=l 

*)^EH)- r(g ;, + 1) 'y . </,<!. (5.4) 

?1=1 

We indeed note that the first term of the above series (5.3)-(5.4) is identical 
to the asymptotic representation of the corresponding functions ^{t) and 
4>{t) of the renewal process of the Mittag-Leffler type, see (4.5)-(4.7). The 



4 The Wright function with parameters A, is defined as 

*\Az) ~Y] rr7T TTT ' x >~ 1 > M gC > zgC - 
— ' nil (An + a) 

n — 

It is an entire function of order p = 1/(1 + A). For detailed information on the functions 
of Wright type the reader may consult e.g. [121 1141 [23l 1241 [27l I49| . and references therein. 
We note that in the classical handbook of the Bateman Project [TD], see Vol. 3, Ch. 
18, presumably for a misprint, the Wright function is considered with A restricted to be 
non-negative. 
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a = 0.25 




p = 0.5 / 






/ p = 0.75 


p = 1 





10 " 
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Figure 2: The functions V(t) (left) and <j>(t) (right) versus t (10~ 2 < t < 10 2 ) 
for the renewal processes of Wright type with j3 = 0.25,0.50,0.75,1. For 
P = 1 the reader would recognize the Box function (extended up to t = 1) 
at left and the delta function (centred in t = 1) at right. 



behaviour near t = of these functions is provided by the first term of their 
asymptotic expansions as t — > 0, namely from 



~ I- A t b/2 exp ( —B r b ^j , <f>(t) ~ C t~ c exp ( -B t 

whera3 



-6 



(5.5) 



2vr(l - 0) /?V(i-/3) 



1/2 



S = (l-/3)/3 b , C 
/? 2-/3 



£1/(1-/3) 

2tt(1 - /3) 



nl/2 



(5.6) 



1-/3' 2(1-/3) 
As far as the functions Vk(t), see (2.14), are concerned, we have 

«o(s) = #(s) = , 



so. 



Vo(t) = ¥(t) 



f 1 (-^) , 0</3<l, 

I e(t) - e(t - 1), p = i, 



(5.7) 



(5. 



5 We take this occasion to point out a misprint (sign error) in the paper |30| . Noting 
that the asymptotic representation as x — + + of the unilateral extremal density of index 
a £ (0, 1) in Eq. (4.15) of [30] is equivalent to ours of cj>(t) as t — > 0, the coefficient ci 
(corresponding to our b in (5-5)- (5-6)), must be taken with the minus sign. 



12 



and 

_ _ e -ks p e -(fc + l)/ 
v k (s) = *(s)[iP(s)] k = — , fc = l,2,..., (5.9) 

from which, in view of the scaling property of the Laplace transform, 

v fc (t) = ( ("^) " (-^) ' 0</?<1 ' (5.10) 

[@{t-k)-@(t-k-l), (3=1. 

Let us close this Section with a discussion on the renewal function m(t). 
Whereas for the process of Mittag-Leffler type we have an explicit expression, 
namely (4.8), we could not find such one for the process of Wright type in 
the case < (3 < 1 from the Laplace transforms (2.19) and (5.2), i.e. 

fh(s) = --, = -s = - Y e ks . 5.11 

s[l-M] s l-e~^ s ti 

In the special case ft = 1 we have (the term by term inversion is allowed!) 

oo 

m(t) = ^e(t-k) = [t], (5.12) 
k=i 

where [t] denotes the greatest integer less than or equal to t. 

In the case < [3 < 1 we do not know how to invert (5.11). We can, 
however, see the asymptotics near zero and near infinity and apply Tauber 
theory: s — > gives fh(s) = l/s 1+/3 , hence, 



rn 



(t) ~ t p /T(l + (3) for t oo , (5.13) 



s —> oo gives m(s) = exp(— s /3 )/s = $(s), hence (if this Tauber trick is 
allowed for such fast decay towards zero), in view of (5.5)-(5.6), 

m(t) ~ ~ ^t fe/2 exp (-Bt 4 ) for t-»0. (5.14) 

6. The Mittag-Leffler and Wright processes in 
comparison 

In this section we intend to compare the renewal processes of Mittag-Leffler 
and Wright type introduced in the Sections 4 and 5, respectively. In this 
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comparison we agree to use the upper indices a and b to distinguish the 
relevant functions characterizing these processes. 

We begin by pointing out the major differences between the survival 
functions \I> a (i) and ^ b (t), provided by Eqs (4.3) and (5.1), respectively, 
for a common index (3 when < (3 < 1. These differences, visible from a 
comparison of the plots (with logarithmic scales) in the left plates of Figures 
1 and 2, can be easily inferred by analytical arguments, as previously pointed 
out (in a preliminary way) in the paper by Mainardi et al. [33j . 

Here we stress again the different behaviour of the two processes in 
the limit (3 — > 1 for t > 0. Whereas ^f a (t) and (f> a (t) tend to the 
exponential exp(— t), ^ b (t) tends to the box function O(i) — Q(t — 1) and 
the corresponding waiting-time pdf 4> b {t) tends to the shifted Dirac delta 
function 5(t — 1). The first of these processes is thus a direct generalization 
of the Poisson process since, for the limiting value (3=1, the Poisson 
process is recovered. In distinct contrast, the second process changes its 
character from stochastic to deterministic: for = 1 at every instant t = n, 
n a natural number, an event happens (and never at other instants) so 
that N{t) = [t], hence trivially m(t) := E(N(t)) = [t] as we had already 
found in (5.12) by summation. We refer to this peculiar counting process 
as the clock process because of its similarity with the tick-tick of a (perfect) 
clock. We also note that in the limit (3 = 1 the densities of the Mittag- 
Leffler and Wright processes have an identical finite first moment since 
p = J °° t exp(-i) dt = / °° t5(t-l)dt = l. 

It is instructive to consider the special value (3 = 1/2 because in cases 
(a) and (b) we have an explicit representation of the corresponding survival 
functions and waiting-time densities in terms of well known functions. 

For the renewal process of Mittag-Leffier type with (3 = 1/2 we have 

y a (t) = E 1/2 (-Vt) = e t exic{Vt) = e l —= \ e~ u du, t>0, (6.1) 

V 71 " J\ft 

where erfc denotes the complementary error function, see e.g. [I], and 

^(t) = -^-E 1/2 (-Vi) = -^=-e t erfc(Vi), t>0. (6.2) 

dt ' yjixt 

For the renewal process of Wright type with (3 = 1/2 we obtain for t > 0, 
*\t) = 1 - $_i /2jl (-lM) = 1 - erfc (J^j = erf (^=) , (6.3) 

and 

4> h {t) = ~ <&-i/ 2 ,o(-lM) = i~ 3/2 exp (-1) . (6.4) 
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We observe that for this particular value of f3 the expression for the densitj{^| 
is obtained not only from the the sum of the convergent series (5.4) but also 
exactly from its asymptotic representation for t — > 0, see (5.5)-(5.6). 

We easily note the common asymptotic (power-law) behaviour of the 
survival and density functions as t — > oo in the cases (a) and (6), that is, 
indicating by the index oo that we mean the leading asymptotic term: 

£-1/2 

Vooit) = —=- , (6.5) 



7T 



t -3/2 

oo (t) = __. (6.6) 

It is now interesting to compare numerically the survival functions and 
the density functions of the two processes, that is (6.1)-(6.2) with (6.3)- 
(6.4), respectively. We find it worthwhile to add in the comparison their 
asymptotic expression (6.5)-(6.6) and also the corresponding functions of 
the Poisson process, namely (3.1)-(3.2). 

In Figure 3, we display the plots concerning the above comparison for 
the survival functions (left plate) and for the density functions (right plate) 
adopting the continuous line for the exact functions of the two long-memory 
processes and the dashed line both for the asymptotic power-law functions 
and for the exponential functions of the Poisson process. The comparison is 
made by using logarithmic scales for 1CP 1 < t < 10 1 , just before the onset 
of the common power law regime; we note that, at least for the case (5 = 1/2 
under consideration, the Wright process reaches the asymptotic power-law 
regime a little bit earlier than the corresponding Mittag-Lemer process. 

In Tables I and II we show for some values of time t the corresponding 
values of the above survival and density functions in comparison, abbrevi- 
ating power law by P-L and Mittag-Lemer by M-L. 



6 We point out C < 4> b (t) = 2^J~ eX ^ \ 4f / ' s J ~ exp(— s 1 ^ 2 ). This Laplace transform 

pair was noted by Levy with respect to the unilateral stable density of order 1/2 and later, 
independently, by Smirnov. In the probability literature the distribution corresponding 
to this pdf is known as the Levy-Smirnov stable distribution. We take this occasion to 
point pout a misprint in the paper [33]. There, in Eq. (3.25) giving the expression of the 
Levy-Smirnov density, the factor 2 in front of was missed. 
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Figure 3: Comparison versus time of the survival functions (left plate) 
and of the corresponding probability densities 4>{t) (right plate). 



Time 


P-L (6.5) 


M-L (6.1) 


Wright (6.3) 


Poisson (3.1) 


0.1 


1.78 10 u 


7.24 lO" 1 


9.74 lO" 1 


9.05 10- 1 


0.5 


7.98 lO" 1 


5.23 lO" 1 


6.83 10" 1 


6.07 10" 1 


1 


5.64 10" 1 


4.28 10" 1 


5.21 10" 1 


3.68 10" 1 


2 


3.99 10" 1 


3.36 10" 1 


3.83 10" 1 


1.35 10" 1 


5 


2.52 10" 1 


2.32 10" 1 


2.48 10" 1 


6.74 10~ 3 


10 


1.78 lO" 1 


1.71 10" 1 


1.77 10" 1 


4.54 10" 5 


20 


1.26 lO" 1 


1.23 10" 1 


1.26 10" 1 


2.06 10- 9 


50 


7.98 10~ 2 


7.90 10" 2 


7.97 10" 2 




100 


5.64 10" 2 


5.61 10- 2 


5.64 10" 2 





Table I: Comparison among the survival functions at different times. 



Time 


P-L (6.6) 


M-L (6.2) 


Wright (6.4) 


Poisson (3.2) 


0.1 


8.92 10 u 


1.06 10- 1 


7.32 lO" 1 


9.05 10- 1 


0.5 


7.98 lO" 1 


2.75 lO" 1 


4.84 10" 1 


6.07 10- 1 


1 


2.82 10" 1 


1.37 lO" 1 


2.20 10" 1 


3.68 10" 1 


2 


9.97 10~ 2 


6.27 10~ 2 


8.80 10" 2 


1.35 10" 1 


5 


2.52 10" 2 


2.00 10" 2 


2.40 10" 2 


6.74 10~ 3 


10 


8.92 10" 3 


7.83 10" 3 


8.70 10" 3 


4.54 10" 5 


20 


3.15 10" 3 


2.94 10~ 3 


3.11 10" 3 


2.06 10" 9 


50 


7.98 10" 4 


7.75 10~ 4 


7.94 10~ 4 




100 


2.82 10" 4 


2.78 10" 4 


2.81 10~ 2 





Table II: Comparison among the density functions at different times. 



We finally consider the functions v^it) for the two long-memory 
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processes. For the functions v%(t) of the Mittag-Lefner process, see (4.10), 
we can take profit of the recurrence relations for repeated integrals of the 
error functions, see e.g. [I], §7.2, pp 299-300, to compute the derivatives of 
the Mittag-Lefner functions in Eqs. (4.10). We recall for n = 0, 1, 2, . . . , 

d 11 / 9 \ 9 

— [e z erfc (z)) = (-l) n 2 n n\ e z I n erfc (z) , (6.7) 

where I n erfc (z) = f™ J n_1 erfc (() d( and J" 1 erfc (z) = 2 exp(-z 2 )/^. 
For the Wright process, see (5.10), we have t>o(i) = ^f b (t) as in (6.3) and 

k = 1,2,... . (6.8) 
7. Renewal processes with reward 

The renewal process can be accompanied by reward that means that at every 
renewal instant a space-like variable makes a random jump from its actual 
position to a new point in "space". "Space" is here meant in a very general 
sense. In the insurance business, e.g., the renewal points are instants where 
the company receives a payment or must give away money to some claim of a 
customer, so space is money. In such process occurring in time and in space, 
also referred to as compound renewal process, the probability distribution of 
jump widths is as relevant as that of waiting times. 

Let us denote by X n the jumps occurring at instants t n , n = 1, 2, 3, . . . , 
and assume that they are iid (real, not necessarily positive) random variables 
with probability density w(x), independent of the waiting time density 4>(t). 
In a physical context the X n s represent the jumps of a diffusing particle (the 
walker), and the resulting random walk model is known as continuous time 
random walk, abbreviated as CTRWS, in that the waiting time is assumed 
to be a random variable with a continuous probability distribution function. 

7 The name CTRW became popular in physics after Montroll, Weiss and Scher (just 
to cite the pioneers) in the 1960's and 1970's published a celebrated series of papers on 
random walks for modelling diffusion processes on lattices, see e.g. [381139 , and the book by 
Weiss [51] with references therein. CTRWs are rather good and general phenomenological 
models for diffusion, including processes of anomalous transport. However, it should 
be noted that the idea of combining a stochastic process for waiting times between two 
consecutive events and another stochastic process which associates a reward or a claim to 
each event dates back at least to the first half of the twentieth century with the so-called 
Cramer-Lundberg model for insurance risk, see [9] for a review. 



v b k (t) 



erfc 



2y/i 



erfc 



f k + 1 
V 2xTt 
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The position x of the walker at time t is 



x(t) 



N(t) 

;(0) + Y. X k 
k=l 



(7.1) 



with N(t) as in (2.5). 

Let us now denote by p(x, t) the probability density (density with respect 
to x) of finding the random walker at the position x at time instant t. We 
assume the initial condition p(x, 0) = 5(x) , meaning that the walker is 
initially at the origin, x(0) = 0. We look for the evolution equation for 
p(x, t) of the compound renewal process. 

Based upon the previous probabilistic arguments we arrive at 



p(x,t) = 5(x)V(t) + / cp(t-t') 
Jo 



+oo 



w{x — x) p{x , t') dx' 



d£ , (7.2) 



called the integral equation of the CTRW. From Eq. (7.2) we recognize the 
role of the survival probability ^>{t) and of the densities cj){t) , w(x) . The 
first term in the RHS of (7.2) expresses the persistence (whose strength 
decreases with increasing time) of the initial position x = 0. The second 
term (a space-time convolution) gives the contribution to p(x, t) from the 
walker sitting in point x' £ R at instant t' < t jumping to point x just at 
instant t , after stopping (or waiting) time t — t' . 

The integral equation (7.2) can be solved by using the machinery of the 
transforms of Laplace and Fourier. Having introduced the notation for the 
Laplace transform in sec. 1, we now quote our notation for the Fourier 
transform, F{f(x); K } = f( K ) = f+™e iKX f(x)dx (k £ R), and for the 
corresponding Fourier convolution of (generalized) functions 



(/i * h) ix) 



+oo 



fi(x') f2(x - x) dx . 



Then, applying the transforms of Fourier and Laplace in succession to 
(7.2) and using the well-known operational rules, we arrive at the famous 
Montroll- Weiss equation, see [39J, 



p(K, 



1 — 4>(s) w(k) 



(7.3) 



As pointed out in [13], this equation can alternatively be derived from 
the Cox formula, see [8] chapter 8 formula (4), describing the process as 
subordination of a random walk to a renewal process. 
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By inverting the transforms one can, in principle, find the evolution 
p(x,t) of the sojourn density for time t running from zero to infinity. In 
fact, recalling that |u5(k)| < 1 and \<j>{s)\ < 1, if k ^ and s / 0, Eq. (7.3) 
becomes 

oo oo 

P(k, s) = §(s) IH S ) w( K )] k = d k(s) w k (K) , (7.4) 

k=0 k=0 

and we promptly obtain the series representation 

oo oo 

p(x, t) = Y Mt) w k {x) = y(t) 5(x) + Y v k(t) w k {x) , (7.5) 

k=0 k=l 

where the functions Vk(t) and wt{x) are obtained by repeated convolutions in 
time and in space, Vk{t) = (^>*(f)* k )(t), see also (2.14), and Wk(x) = (w* k )(x), 
respectively. In particular, v (t) = (^ * 5)(t) = v\(t) = (* * 4>){t), 

wo(x) = 5(x), w\(x) = w(x). In the R.H.S of Eq (7.5) we have isolated the 
first singular term related to the initial condition p(x, 0) = ^(0) 5(x) = 5(x). 

A special case of the integral equation (7.2) is obtained for the compound 
Poisson process where <j)(t) = e~ l (as in (3.2), with A = 1 for simplicity). 
Then, the corresponding equation reduces after some manipulations, that 
best are carried out in the Laplace-Fourier domain, to the Kolmogorov- Feller 
equation: 

— p(x,t) = —p(x,t)+ / w(x — x')p(x , t) dx' , (7-6) 
at J-oo 

which is the master equation of the compound Poisson process. In view of 
Eqs (3.5) and (7.5) the solution reads 

oo 

p(x,t) = Y,Ti e ~ t M x )- ( 7 - 7 ) 

k=0 

When the survival probability is the Mittag-Leffler function introduced in 
(4.3), the master equation of the corresponding compound process of Mittag- 
Leffler type can be shown to be 

/+oo 
w(x-x')p(x',t)dx , 0</3<l, (7.8) 
-oo 

where tD* denotes the time fractional derivative of order [3 in the Caputo 
sense. For a (detailed) derivation of Eq (7.8 ) we refer to the paper by 
Mainardi et al. [33], in which the results have been obtained by an approach 
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different from that adopted in a previous paper by Hilfer and Anton [22]. In 
this case, in view of Eqs (4.10) and (7.5), the solution of the master equation 
(7.8) reads: 



p(x, t) = E 



k=l 



Ef\-t p )w k {x), </?<!. (7.9) 



When the survival probability is the Wright function introduced in (5.1), 
we resist the temptation to work out an explicit form for the the master 
equation of the corresponding compound process of Wright type, but we 
content ourselves to give its series solution derived from (7.5) using Eqs 
(5.8), (5.10), as 



p(x,t) 



1 - 



1 

~tP 



S(x) + 



(7.10) 



oo 

E 

fc=l 



13,1 



k 



fe + 1 
~t?~ 



w k (x), <(3<1. 



Recalling what we have observed in Section 6 on the contrasting behaviours 
of the Mittag-Lefher process and the Wright process in the limit (3 = 1 
we recognize that with rewards the first of these processes becomes the 
compound Poisson process ((7.9) with (3 = 1 reduces to (7.7)), whereas the 
second goes over into the process, whose sojourn pdf is directly obtained 
from the Laplace inversion of Eqs (5.7) and (5.9) for (3 = 1. Recalling the 
natural convention Wq(x) = 5(x), we obtain 



p(x, t) = J2 -k)-e(t-k- 1)] w k (x) . 

k=0 



(7.11) 



This simply means: p(x,t) = w k (x) for k < t < k + 1 (k = 0,1,2,...), 
or, more concisely, p(x,t) = wm(x). This limiting process is nothing but a 
simple random walk, discrete in time and Markovian on the set of integer 
time-values: at every instant [t] a jump occurs with pdf w(x), hence [t] 
jumps have occurred exact up to instant t, which leads to [t]-fold convolution 
W[t\(x). We have rigorously investigated such random walks in [T7j. However, 
in the framework of the CTRW theory (where the resulting random walk, 
discontinuous or continuous in space, is implicitly intended to be continuous 
in time), this case appears a peculiar (non-Markovian) process, as discussed 
by Weiss [51 j in his book at p. 47, see also [33] . 
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8. Time fractional diffusion as limit of CTRW 

There are several ways to generalize the classical diffusion equation by 
introducing space and/or time derivatives of fractional order. We mention 
here the seminal papers by Schneider & Wyss [16] of 1989 for the time 
fractional diffusion equation and the likewise influential paper by Saichev & 
Zaslavsky [42J of 1997 for diffusion fractional in time as well in space. In the 
recent literature several authors stress the viewpoint of subordination, see 
e.g. [21 [32l [351 [52], an d special attention is being paid to diffusion equations 
with distributed orders of fractional temporal or/and spatial differentiation, 
see e.g. [6l [191 EH]- The transition from CTRW to such generalized 
types of diffusion has been investigated by different methods, not only from 
mathematical interest but also due to applications in Physics, Chemistry (for 
generalized Fokker-Planck and Liouville equations see e.g. [3l 14"! |2~T | 122" 1 [47] ) . 
and other Applied Sciences including Economics and Finance (for financial 
markets see e.g. [20l [33l Ej). For a good list of references on these topics 
we recommend the review papers by Metzler and Klafter [361 137j . 

Gorenflo and Mainardi, see e.g. [13 [19], have, under proper conditions 
on the tails of the probability distributions for jumps and waiting times, 
investigated for Eq. (7.2) the so-called well scaled transition to the diffusion 
limit. This transition is obtained by making smaller all jumps by a positive 
factor h and all waiting times by a positive factor r related to h by a proper 
scaling relation, and then letting h and r tend to zero. An alternative 
interpretation is that we look at the same process from far away and after 
long time, so that spatial distances and time intervals of normal size appear 
very small but always related through proper scaling. In this limit the 
integral equation (7.2) has been shown to reduce to a partial differential 
equation with fractional derivatives in space and/or in time, referred to as a 
space-time fractional diffusion equation. This is also the topic of the recent 
paper by Scalas et al. [15] and, in a more general framework, of the paper 
by Gorenflo and Abdel-Rehim [13]. For a variant of well-scaled transition 
(there not called so) see [50]. In this paper we limit ourselves to recall the 
well scaled transition from the CTRW integral equation (7.2) to the time 
fractional diffusion equation (TFDE), 



with < P < 1, subjected to the initial condition u(0,t) = 5(x). For 
(3 = 1 we recover the standard diffusion equation, for which the fundamental 
solution is the Gaussian probability density evolving in time with variance 
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a 2 = 2t, 

n(x , t) = J_ t -i/2 e -x 2 /(4t). (8 . 2) 

2y/7T 

This probability law is known to govern the classical Brownian model for 
the phenomenon of normal diffusion. Also for < (3 < 1 the fundamental 
solution of t (8.1) can be interpreted as a (spatially-symmetric) density 
evolving in time, but exhibiting stretched exponential tails with a variance 
proportional to t@ , implying a phenomenon of slow anomalous diffusion. 
We have, see for details pU CDS EU E3 EUJ M\ , 

u( x ,t) = ^M m (\x\/t /3 / 2 ), (8.3) 

where 

^' s S ^g-w <8 ' 4) 



•7' 



\/t^/ 2 acting as a similarity variable. We stress that, according to footnote 
( 4 ), Mp/ 2 {y) = $-/3/24-j8/2(z/)) hence the fundamental solution of the TFDE 
is a Wright function! The TFDE (8.1) can be derived from the CTRW 
integral equation (7.2), by properly rescaling the waiting time and the jump 
widths and passing to the diffusion limit. Making smaller all waiting times 
by a positive factor r , all jumps by a positive factor h , we get, for n S N, the 
jump instants t n {r) = tT\ + tT 2 + . . . + rT n and the jump sums, xo(h) = 0, 
x n (h) = hXi + hX 2 + . • • + hX n . The reduced waiting times rT n all have 
the pdf (j) T (t) = 0(t/r)/r, t > 0, analogously the reduced jumps hX n all 
have the pdf Wh(x) = w{x/h)/h, x £ R. The probability density p(x,t) so 
obtained is denoted by p^ T (x,£). Then, the transition to the diffusion limit 
consists in sending r — ► and h — * under an appropriate relation between 
r and h and deriving the evolution equation satisfied by the limiting pdf 



8 In his first 1993 analysis of the time fractional diffusion equation, one of the present 
authors (F.M.) [25] introduced the two auxiliary functions (of Wright type), F v (z) := 
< l > -y,o(— z) and M v {z) := $—v,i—v(—z) with < v < 1, inter-related through F v {z) = 
v z M v (z) . Being in that time only aware of the Bateman project where the parameter 
A of the Wright function &\ tll (z) was erroneously restricted to non-negative values, F.M. 
thought to have extended the original Wright function. It was just Professor Stankovic 
during the presentation of the paper [31] in the Conference Transform Methods and Special 
Functions, Sofia 1994, wno informed F.M. that this extension for —1 < /i < was already 
made just by Wright himself in 1940 (following his previous papers in 1930's). On this 
special occasion F.M. wants to renew his personal gratitude to Professor Stankovic for 
this earlier information, that has induced him to study the original papers by Wright and 
work (also in collaboration) on the functions of the Wright type for further applications. 



22 



Let us now resume the approach by Gorenflo and Mainardi, for giving 
the conditions on w(x), (j>(t) and the scaling relation between h and r, which 
allow po o(x,t) to be identified with the fundamental solution u(x,t) of our 
TFDE (8.1). Assuming the jump pdf w(x) to be symmetric with finite 
variance a 2 and the waiting time pdf <fi(t) either with finite mean p (relevant 
in the case (3 = 1), or, with c > and some (3 € (0, 1) , <p{t) ~ ct~^ +1 ) for 
t -> oo, and setting p = a 2 , A = p if (3 = 1, A = cT(l - (3)/ (5 if < (3 < 1, 
the required scaling relation for the diffusion limit is At' 3 = pti 2 . 

We now will compare, at fixed times, the spatial probability distributions 
provided by the fundamental solutions of the integral equation of the 
CTRW as given by the series (7.5) with the fundamental solutions of the 
corresponding limiting TFDE as given by the function (8.3) of Wright type. 
This will enable us, for the compound Poisson, Mittag-Leffier and Wright 
processes, to investigate from a numerical view the increasing quality of 
approximation to the diffusion limit with advancing time, independently 
from the scaling relation. In this numerical comparison, in order to avoid the 
singular terms originated by delta functions, we prefer to consider, instead of 
the spatial densities p(x,t), w(x) and u(x,t), the corresponding probability 
cumulative functions that we denote by capital letters, namely P(x, t), W(x) 
and U(x, t). Thus, by integrating in space (7.5) from — oo to x, and denoting 
by Q(x) the unit step Heaviside function, we have 



p(x,t) = tf(t)e(z) + X> fc (t)w- fc (z) 
k=i 



(8.5) 



In an analogue way, for the TFDE we have 



U(x,t) = - 



rx/tl 3 / 2 

l + y o M p/2 {y)dy 



O<0<1, 



(8.6) 



that in the case (3=1 reduces to 



U{x,t) = \ 



1 + erf 



x 



2Vi 



(8.7) 



In our case studies we have chosen t = 1 and t = 10, as typical short 
and long times. As waiting time distributions we have used the exponential 
distribution for the compound Poisson process so that the v^it) are given 
by (3.5), and the Mittag-Leffier and Wright distributions for (3 = 1/2, for 
the corresponding compound processes so that v^t) are given by (4.10) with 
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(6.7), and (6.8), respectively. For all the processes we have used a symmetric 
Gaussian for the common jump probability distribution, for which we have 



1 e -*/( 4fc ) 1 

Wk{x) = V? Vk ' Wk{x) = 2 



1 + erf 



In figures 4,5,6, regarding respectively, the compound Poisson, Mittag- 
Leffler and Wright processes, we report in continuous line the pcf for the 
CTRW whereas in dashed line the corresponding cumulative functions of 
the limiting TFDE. The values of the survival function ^(t) for t = 1 and 
t = 10 are given in Table I; however, they are clearly visible as the height of 
the vertical part (at x = 0) of the continuous line. 

We note that for j3 = 1/2 there is no substantial difference between 
the cases of Mittag-Leffler and Wright compound processes: at these times 
the dominant effect is due to the common asymptotic behaviour of the 
corresponding distributions resulting from the diffusion limit. In other 
words: the results are consistent with the asymptotic properties of the 
diffusion limit. We also note that in this comparison we cannot take direct 
profit of the scaling relation of the diffusion limit, because the fundamental 
solution of the CTRW integral equation has no self-similarity property at 
variance to the fundamental solution of the corresponding TFDE. 

We point out that also Barkai [3] has investigated the diffusion limit for 
the CTRW process of Wright type but considering different case studies and 
using different notation and terminology. 



9. Conclusions 

After sketching the basic principles of renewal theory we have discussed the 
classical Poisson process. It is well known that in this process the sum of 
k waiting times obeys the so-called Erlang distribution, named after the 
Danish telecommunication engineer Agner Krarup Erlang (1878-1929 )q 

Because the exponential function is just a special case of the Mittag- 
Leffler function, one may expect that an analogous explicit representation 
exists for the situation of waiting times distributed according to the Mittag- 
Leffler law Vl/(t) = Ep^-t 13 ), < (3 < 1, with denoting the survival 
function (the probability of no event occurring in an interval of length t just 
after an event). We have thus presented such explicit representations, for 
the above renewal process that we have called renewal process of Mittag- 
Leffler type. We have also considered another process where the survival 

9 See http://plus.maths.org/issue2/erlang/ for information on Erlang. 
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function, being related to the unilateral Levy stable law of index (3, is given 
as an expression involving the Wright function. This process has been called 
renewal process of Wright type. 

In both of these renewal processes the waiting time probability 
distribution is long-tailed, and there is a parameter (3, < (3 < 1. The first 
of these processes is a direct generalization of the Poisson process as for the 
limiting value (3 = 1 the Poisson process is recovered. In distinct contrast, 
the second process changes its character from stochastic to deterministic: 
for j3 = 1 at every instant t = n, n a natural number, an event happens 
(and never for other instants). Then, we have reviewed the general theory 
of renewal processes with reward, the so called compound renewal processes, 
known in physics and chemistry literature as continuous time random walks. 

Using the fact that the sojourn probability density of a compound 
renewal process can be written as an infinite series involving the convolution 
powers of the waiting times and the rewards (the jumps), we have simulated 
these processes (all three cases: first Poisson, then, with (3 = 1/2, Mittag- 
Leffter and Wright). This comparison has allowed us to make visible how in 
the large time regime the behaviour of the probability cumulative functions 
approaches that of the corresponding functions for time-fractional diffusion 
processes, to which our compound renewal processes are known to weakly 
converge in the well scaled transition to the diffusion limit. Nicely one can 
see that the survival function decays towards zero very fast (exponentially) in 
the Poisson process, but slowly (like a power with negative exponent) in the 
other two processes. Incidentally we have also stressed, from a theoretical 
view-point, the remarkable contrast in the behaviour of the non-Markovian 
compound processes of Mittag-Leffler and of Wright type in the limit (3 = 1, 
yielding the classical compound Poisson process (which is continuous in time 
and Markovian) for the first type, but a process discrete in time (Markovian 
on the set of the integer time-values, but non Markovian on the set of the 
real time-values) for the second. 
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